source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/SLE_TRPM2_MfMo/'
kn.pal2 <- c('#d44','#485ffc')

# washburn 25 ---------
washburn <- read_rds('mission/SLE_TRPM2_MfMo/GSE290695_ileum_treatment_naive_CD.rds')

washburn |> DimPlot()

washburn |> DotPlot('TRPM2', assay = 'RNA', cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Cell type',
       title = "TRPM2 expression in Crohn's disease ileum tissue",
       subtitle = 'GSE290695 (n=27)')

washburn$cell_typev2 

sobj_mfmo <- washburn |>
  filter(str_detect(cell_typev2, 'Monocyte|Macrop'))

sobj_mfmo <- sobj_mfmo |> FindClusters()

sobj_mfmo |> DotPlot('TRPM2', assay = 'RNA')

sobj_mfmo |>
  filter(str_detect(cell_typev2, 'Macro')) |>
  DotPlot('TRPM2', assay = 'RNA', cols = 'RdYlBu', cluster.idents = T)

last_plot() |>
  pluck('data') |>
  BubblePlot(size = 5) +
  labs(x = 'Gene', y = 'Cluster',
       title = "TRPM2 expression in Crohn's disease ileum tissue macrophage",
       subtitle = 'GSE290695 (n=27)') 

## corr with TRPM2 ---------
plsm_membr_go <- map_go_gene('GO:0005886')

ibd_louv08_mean <- sobj_mfmo |>
  AverageExpression(features = plsm_membr_go$SYMBOL, assays = 'RNA') |>
  pluck('RNA') |>
  t() |>
  as_tibble(rownames = 'cluster')

wb_louv08_m2 <- ibd_louv08_mean |>
  dplyr::select(cluster, TRPM2)

ibd_m2_cor1 <- ibd_louv08_mean |>
  pivot_longer(-1, names_to = 'gene') |>
  left_join(wb_louv08_m2) |>
  summarize(correlation = cor(value, TRPM2),
            p.val = cor.test(value, TRPM2)$p.value,
            .by = gene) |>
  na.omit()

ibd_m2_cor1 |>
  mutate(gene = fct_reorder(gene, correlation)) |>
  slice_max(gene, n = 20) |>
  ggplot(aes(correlation, gene, fill = p.val)) +
  geom_point(shape = 21) +
  scale_fill_distiller(palette = 'Reds') +
  theme_bw(base_size = 11, base_family = 'ArialMT') +
  labs(title = 'Top 20 membrane protein gene correlated with TRPM2 in CD monocyte',
       fill = 'p value') +
  theme_jpub

publish_source_plot('top20.membrane.gene.cor.TRPM2.mono', width = 60)

sle_m2_cor <-
read_csv('mission/SLE_TRPM2_MfMo/results/SLE.PBMC.membr.gene.cor.TRPM2.mono.csv')

ibd_sle_m2_cor <- ibd_m2_cor1 |>
  inner_join(sle_m2_cor, join_by(gene))

ibd_sle_head <- ibd_sle_m2_cor |>
  mutate(gene = fct_reorder(gene, correlation.x + correlation.y)) |>
  slice_max(gene, n = 6)
  
ibd_sle_m2_cor |>
  ggplot(aes(correlation.x, correlation.y)) +
  geom_point(alpha = .1) +
  geom_point(data = ibd_sle_head) +
  geom_text_repel(data = ibd_sle_head, aes(label = gene)) +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  stat_cor() +
  theme_bw() +
  stat_smooth(method = 'lm') +
  labs(title = 'Membrane protein genes expression correlation with TRPM2 in monocytes',
       x = 'Correlation in Crohn disease', y = 'Correlation in SLE')

# Tabula muris M2 cor ------
plsm_membr_marker <- read_csv('mission/SLE_TRPM2_MfMo/results/plsm_membr_marker.csv')

## bone marrow -----------
mm_marrow <- load('~/droplet_Marrow_seurat_tiss.Robj')
mm_marrow
mm_marrow <- tiss |> UpdateSeuratObject()
mm_marrow |> VlnPlot('percent.ribo', pt.size = 0)
mm_marrow |> DimPlot(group.by = 'cell_ontology_class')

Idents(mm_marrow) <- mm_marrow$cell_ontology_class

mm_marrow |> DotPlot('Trpm2', cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Trpm2 expression in healthy mouse bone marrow',
       subtitle = 'Tabula Muris, Nature, 2018') +
  theme_jpub

bm_macro <- mm_marrow |>
  filter(str_detect(cell_ontology_class, 'macro'))

bm_macro <- bm_macro |>
  quick_process_seurat(batch = 'channel', skip_norm = T)

bm_macro |> DotPlot('Trpm2', scale = F, cols = 'RdYlBu')

bm_macro |> FindMarkers(ident.1 = 1, features = 'Trpm2')

bm_macro_pm <- bm_macro |>
  get_abundance_sc_long(features = plsm_membr_marker$SYMBOL)

bm_macro_pm <- bm_macro_pm |>
  summarise(sum = sum(.abundance_RNA), .by = .feature) |>
  filter(sum > 0) |>
  mutate(sum = NULL) |>
  left_join(bm_macro_pm)

bm_macro_m2 <- bm_macro_pm |>
  filter(.feature == 'Trpm2') |>
  pivot_wider(names_from = .feature, values_from = .abundance_RNA)

bm_macro_cor <- bm_macro_pm |>
  left_join(bm_macro_m2) |>
  summarize(correlation = cor(.abundance_RNA, Trpm2),
          p.val = cor.test(.abundance_RNA, Trpm2)$p.value,
          .by = .feature) |>
  na.omit() |>
  mutate(rank = rank(correlation)) |>
  arrange(rank)

bm_macro_cor_min <- bm_macro_cor |>
  slice_min(correlation, n = 5)

bm_macro_cor_mm <- bm_macro_cor |>
  slice_max(correlation, n = 6) |>
  bind_rows(bm_macro_cor_min)

bm_macro_cor |>
  mutate(.feature = fct_reorder(.feature, correlation)) |>
  ggplot(aes(correlation, rank)) +
  geom_point(size = .01) +
  geom_point(data = bm_macro_cor_mm, size = 1, color = 'red') +
  geom_label_repel(data = bm_macro_cor_mm, aes(label = .feature), size = 2,
                   nudge_y = ifelse(bm_macro_cor_mm$correlation > 0, -300, 300)) +
  theme_bw(base_size = 6) +
  labs(title = 'Membrane protein gene correlated with TRPM2 in mouse BM',
       fill = 'p value') +
  theme(plot.title.position = 'plot', legend.key.size = unit(4, 'mm'))

publish_source_plot('mouse.BM.membr.gene.cor.TRPM2.mono', width = 70)

## spleen --------
mm_spleen <- load('~/droplet_Spleen_seurat_tiss.Robj')

mm_spleen <- tiss |> UpdateSeuratObject()

mm_spleen |> DimPlot(group.by = 'cell_ontology_class')

Idents(mm_spleen) <- mm_spleen$cell_ontology_class

mm_spleen |> DotPlot('Trpm2', cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Trpm2 expression in healthy mouse spleen',
       subtitle = 'Tabula Muris, Nature, 2018') +
  theme_jpub

sp_macro <- mm_spleen |>
  filter(str_detect(cell_ontology_class, 'macro'))

sp_macro <- sp_macro |>
  quick_process_seurat(batch = 'channel', skip_norm = T)

sp_macro |> DotPlot('Trpm2', cols = 'RdYlBu') +
  labs(x = 'gene', y = 'cluster',
       title = 'Trpm2 expression in healthy mouse spleen macrophage')

sp_leiden_mean <- sp_macro |>
  AverageExpression(features = str_to_title(plsm_membr_marker$SYMBOL))

sp_leiden_mean <- sp_leiden_mean$RNA |>
  t() |>
  as_tibble(rownames = 'cluster')

sp_leiden_m2 <- sp_leiden_mean |>
  dplyr::select(cluster, Trpm2)

sp_m2_cor <- sp_leiden_mean |>
  pivot_longer(-1, names_to = 'gene') |>
  left_join(sp_leiden_m2) |>
  summarize(correlation = cor(value, Trpm2),
            p.val = cor.test(value, Trpm2)$p.value,
            .by = gene) |>
  na.omit() |>
  mutate(rank = rank(correlation)) |>
  arrange(rank)

sp_macro_cor_min <- sp_m2_cor |>
  slice_min(correlation, n = 5)

sp_macro_cor_mm <- sp_m2_cor |>
  slice_max(correlation, n = 6) |>
  bind_rows(sp_macro_cor_min)

sp_m2_cor |>
  mutate(gene = fct_reorder(gene, correlation)) |>
  ggplot(aes(correlation, rank)) +
  geom_point(size = .01) +
  geom_point(data = sp_macro_cor_mm, size = 1, color = 'red') +
  geom_label_repel(data = sp_macro_cor_mm, aes(label = gene), size = 2,
                   nudge_y = ifelse(sp_macro_cor_mm$correlation > 0, -300, 300)) +
  theme_bw(base_size = 6) +
  labs(title = 'Membrane protein gene correlated with TRPM2 in mouse spleen',
       fill = 'p value') +
  theme(plot.title.position = 'plot', legend.key.size = unit(4, 'mm'))

publish_source_plot('mouse.spleen.membr.gene.cor.TRPM2.macro', width = 70)

# GSE266616 ----------
meta <- GEOquery::getGEO('GSE266616')

meta$`GSE266616-GPL18573_series_matrix.txt.gz`

meta81 <-
meta$`GSE266616-GPL24676_series_matrix.txt.gz` |>
  pData() |>
  as_tibble(.name_repair = 'universal')

group81 <- meta81 |>
  select(title, condition.ch1)

ibd_path <-
list.dirs('mission/SLE_TRPM2_MfMo/zhao25ibd/')

ibd_paths <-
tibble(path = ibd_path,
       name = str_remove(path, '.+\\/') |> str_remove('_corrected') |>
         str_replace('-', '_')) |>
  inner_join(group81, join_by(name == title))

tibble(path = ibd_path,
       name = str_remove(path, '.+\\/') |> str_remove('_corrected') |>
         str_replace('-', '_')) |>
  full_join(group81, join_by(name == title)) |>
  DT::datatable()

ibd_mex <-
ibd_paths |>
  pull(path, name = name) |>
  Read10X()

ibd_mex |> dim()

# 484k cells
sobj <- CreateSeuratObject(ibd_mex, min.cells = 3, min.features = 200)

sobj <- sobj |>
  PercentageFeatureSet('^MT-', col.name = 'mito_ratio')

sobj |> VlnPlot('mito_ratio', pt.size = 0)

# 175k cells
sobj <- sobj |> filter(mito_ratio < 10)

sobj <- sobj |> mutate(orig.ident = str_extract(.cell, '.+(?=_)'))

Idents(sobj) <- sobj$orig.ident

sobj <- group81 |>
  mutate(orig.ident = title, group = condition.ch1, .keep = 'none') |>
  left_join(x = sobj, y = _)

sobj <- sobj |>
  quick_process_seurat()

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(hpca, new_label = 'hpca_main')

sobj |>
  DimPlot(group.by = 'hpca_main', cols = 'Paired')

sobj |>
  DotPlot(c('C1QA','CD14', 'TRPM2'), cols = 'RdYlBu', cluster.idents = T)

FindMarkers(sobj, ident.1 = 13, only.pos = T) |>
  as_tibble(rownames = 'gene')

sobj |>
  DotPlot(cc.genes.updated.2019)

eec_mrk <- c('CHGA','SYP','SCG2','NEUROD1','PAX6')
paneth <- c('DEFA5','DEFA6','MMP7','PLA2G2A','SOX9')
goblet <- c('MUC2','TFF3','CLCA1')

cell_markers <- list(
  # Epithelial Cells
  Enterocyte = c("FABP1", "SI", "APOA4"),
  Paneth = c("LYZ", "DEFA5", "DEFA6"),
  Goblet = c("MUC2", "TFF3", "FCGBP"),
  Enteroendocrine = c("CHGA", "GCG", "SST"),
  Tuft = c("DCLK1", "TRPM5", "PTGS1"),
  Stem = c("LGR5", "ASCL2", "OLFM4"),
  
  # Immune Cells
  B_cell = c("CD19", "MS4A1", "CD79A"),
  T_cell = c("CD3D", "CD3E", "CD8A"),
  Macrophage = c("CD68", "CD163", "CSF1R"),
  Dendritic = c("CD1C", "CLEC9A", "CD83"),
  Neutrophil = c("S100A8", "S100A9", "FCGR3B"),
  
  # Stromal Cells
  Fibroblast = c("COL1A1", "COL3A1", "PDGFRA"),
  Myofibroblast = c("ACTA2", "TAGLN", "MYH11"),
  Endothelial = c("PECAM1", "VWF", "CDH5"),
  Pericyte = c("PDGFRB", "RGS5", "CSPG4")
)

sobj |> DotPlot(cell_markers, cols = 'RdYlBu', cluster.idents = T) +
  RotatedAxis()

sobj <- sobj |>
  mutate(cell_type = case_match(as.numeric(seurat_clusters),
                                c(11,17) ~ 'Epithelial_cells',
                                29 ~ 'Neutrophils',
                                25 ~ 'EEC',
                                20 ~ 'Tuft_cells',
                                8 ~ 'Macrophages',
                                13 ~ 'Mast_cells',
                                c(1,4,5,19) ~ 'Plasma_cells',
                                7 ~ 'T_cells',
                                .default = hpca_main))

paired.pal <- RColorBrewer::brewer.pal(name = 'Paired', n = 12)

paired.pal <- c(paired.pal, '#AEC7E8', '#FFD7BE', '#47627A')

sobj |>
  DimPlot(group.by = 'cell_type', cols = paired.pal, raster = T) +
  ggtitle('PBMC') +
  theme_jpub +
  theme_classic(base_size = 6, base_family = 'ArialMT')

publish_pdf('IBD.PBMC.umap.pdf', width = 80)

sobj <- sobj |>
  mutate(group = ifelse(group == 'CD', 'IBD', 'HC') |>
           fct_relevel('HC'))

## rds checkpoint ----------
sobj |> write_rds('mission/SLE_TRPM2_MfMo/zhang_ibd.rds')

sobj <- read_rds('mission/SLE_TRPM2_MfMo/zhang_ibd.rds')

## TRPM2 bubbleplot --------
m2_2d <- sobj |>
  DotPlot2d('TRPM2',group.x = group, group.y = cell_type) |>
  pluck('data')

m2_2d |>
  mutate(group.x = ifelse(group.x == 'CD', 'IBD', 'HC'),
         group.y = fct_reorder(group.y, avg.exp)) |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type',
       title = 'TRPM2 expression in colon tissue') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('IBD.PBMC.TRPM2.dotplot')

## TRPM2 featureplot --------
sobj |> FeaturePlot('TRPM2', cols = c('lightgrey','red'), order = T,
                    split.by = 'group') &
  theme_jpub

## TRPM2 mRNA IBD vs HC logfc --------
m2_ibdvhc <-
sobj |> FindMarkersAcrossVar(split.by = 'cell_type', group.by = 'group',
                             ident.1 = 'IBD', features = 'TRPM2')

sobj_mf <-
sobj |> filter(cell_type == 'Macrophages')

## TRPM2 in MF by sample -------
sobj_mf |> DotPlot('TRPM2', cols = 'RdYlBu', group.by = 'orig.ident')

m2_by_sample <- last_plot() |>
  pluck('data')

m2_by_sample |>
  left_join(ibd_paths, join_by(id == name)) |>
  filter(avg.exp > 2) |>
  mutate(tissue = str_extract(id, 'AC|TI')) |>
  ggplot(aes(condition.ch1, avg.exp)) +
  geom_boxplot() +
  geom_jitter(width = .1, height = 0) +
  theme_pubr() +
  facet_wrap(~tissue) +
  stat_compare_means(method = 't.test')

sobj_mf <- sobj_mf |>
  quick_process_seurat(skip_norm = T)

myeloid_markers <- list(
  Macrophage = c("CD68", "CD163", "MSR1"),
  Monocyte = c("CD14", "FCGR3A", "S100A8", "C1QA"),
  Dendritic = c("CD1C", "CLEC9A", "CD83"),
  Neutrophil = c("CSF3R", "FCGR3B", "S100A9"),
  pDC = c('LILRA4','CLEC4C','TCF4'),
  'TRPM2'
)

sobj_mf |>
  DotPlot(myeloid_markers, cols = 'RdYlBu', cluster.idents = T) +
  RotatedAxis()

sobj_myl <- sobj_mf |>
  mutate(fine_type = case_match(as.numeric(seurat_clusters),
                                19 ~ 'pDC',
                                c(4,17,14,15,20) ~ 'cDC',
                                c(9,3) ~ 'Neutrophil',
                                .default = 'Macrophage'))

sobj_myl |> DimPlot(group.by = 'fine_type', cols = 'RdYlBu')

sobj_mf <- sobj_myl |>
  filter(fine_type == 'Macrophage')

sobj_mf <- sobj_mf |>
  mutate(group = ifelse(group == 'CD', 'IBD', 'HC'))

## MF rds ---------
sobj_mf |> write_rds('mission/SLE_TRPM2_MfMo/zhang_ibd_MF.rds')

sobj_mf <- read_rds('mission/SLE_TRPM2_MfMo/zhang_ibd_MF.rds')

## TRPM2 in MF by sample -------
sobj_mf |> DotPlot('TRPM2', cols = 'RdYlBu', group.by = 'orig.ident')

m2_by_sample <- last_plot() |>
  pluck('data')

m2_by_sample |>
  left_join(ibd_paths, join_by(id == name)) |>
  mutate(group = ifelse(condition.ch1 == 'CD', 'IBD', 'HC')) |>
  #filter(avg.exp < 2) |>
  #mutate(tissue = str_extract(id, 'AC|TI')) |>
  ggplot(aes(group, avg.exp, color = group)) +
  geom_boxplot(outliers = F) +
  geom_jitter(width = .1, height = 0) +
  theme_pubr() +
  #facet_wrap(~tissue) +
  stat_compare_means(method = 't.test', comparisons = list(c('HC','IBD')),
                     color = 'black', label = 'p.signif') +
  labs(x = 'Group', y = 'Average expression',
       title = 'TRPM2 expression in colon tissue macrophages')

## MF M2-dotplot & UMAP ---------
sobj_mf <- sobj_mf |> quick_process_seurat(skip_norm = T)

sobj_mf <- sobj_mf |> FindClusters(algorithm = 4, resolution = .4,
                                   random.seed = 1)

sobj_mf |> DimPlot(cols = paired.pal, label = T, label.box = T, repel = T)

sobj_mf |>
  DotPlot2d('TRPM2', seurat_clusters, group)

sobj_mf |>
  filter(group == 'IBD') |>
  DotPlot('TRPM2', cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Macrophage clusters')

mf_leiden_m2 <- last_plot() |>
  pluck('data')

mf_leiden_m2 |>
  mutate(seurat_clusters = id, avg.exp.scaled, avg.exp, .keep = 'none') |>
  left_join(x = sobj_mf, y = _) |>
  ggplot(aes(umap_1, umap_2, color = avg.exp.scaled)) +
  geom_point(size = 1) +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_classic() +
  labs(color = 'Mean expr',
       title = 'Mean expression of TRPM2 in macrophages clusters')

AutoPointSize(tibble(x = rnorm(n=10300)))

## define M2-hi MF -------
sobj_mf <- sobj_mf |>
  mutate(trpm2_type = ifelse(seurat_clusters %in% c(2,4,8,10),
                             'TRPM2-hi MF', 'TRPM2-lo MF'))

sobj_mf |> DimPlot(group.by = 'trpm2_type')

## inflam score --------
flare.gsebp.tb <-
  read_csv('mission/SLE_TRPM2_MfMo/results/flare.m2hvl.gsebp.csv')

inflam.gene <- flare.gsebp.tb |>
  filter(str_detect(Description, '^inflamma')) |>
  pull(core_enrichment) |>
  str_split_1('/')

sobj_mf <- sobj_mf |>
  AddModuleScore(features = list(inflam.gene), name = 'inflam')

inflam.module <- sobj_mf |>
  filter(group == 'IBD') |>
  summarise(inflam.score = mean(inflam1), .by = seurat_clusters) |>
  mutate(id = seurat_clusters, .keep = 'unused')

mf_leiden_m2 |>
  mutate(subtype = ifelse(id %in% c(2,4,8,10), 'TRPM2-hi MF', 'TRPM2-lo MF')) |>
  left_join(inflam.module) |>
  filter(id != 8) |>
  ggplot(aes(avg.exp, inflam.score)) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2, label.x = .2, label.y = .2) +
  labs(x = 'Average expression of TRPM2', y = 'Inflammatory signature score',
       title = 'IBD colon macrophage clusters') +
  theme_jpub

publish_source_plot('IBD.colon.MF.inflam.score-M2.cor.no8', width = 70)

## M2-hi/lo MF frac ----------
mf_frac_by_sample <- sobj_mf |>
  filter(seurat_clusters != 8) |>
  calc_frac_conf_on_grouped_count(orig.ident, trpm2_type)

sobj_mf |>
  distinct(orig.ident, group) |>
  left_join(mf_frac_by_sample) |>
  ggplot(aes(group, fraction, color = group)) +
  geom_boxplot() +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~trpm2_type) +
  theme_bw() +
  scale_color_manual(values = kn.pal2) +
  stat_compare_means(comparisons = list(c('HC', 'IBD')), method = 't.test') +
  labs(title = 'Macrophage subsets proportion in IBD colon tissue')

## inflam violin ---------
top15.inflam <-
  read.csv('mission/SLE_TRPM2_MfMo/results/top15.inflamm.gene.csv')

sobj_mf |>
  filter(group == 'IBD', seurat_clusters != 8) |>
  bill.violin(top15.inflam$gene, group.by = trpm2_type, facet.ncol = 5) +
  labs(x = 'Cell type', y = 'Normalized expression', fill = 'Cell type',
       title = 'Inflammatory response genes in IBD colon') +
  scale_fill_hue(labels = c('TRPM2-hi Macro', 'TRPM2-lo Macro')) +
  theme_classic() +
  theme_jpub +
  theme(axis.text.x = element_blank(), legend.position = 'top')

## M2-hi GSEA -------
library(clusterProfiler)
m2hvl_deg <- sobj_mf |>
  FindMarkersAcrossVar(split.by = 'group', group.by = 'trpm2_type',
                       ident.1 = 'TRPM2-hi MF')

m2hvl_deg |> write_source_csv('ibd.colon.mf.m2hvl.deg')

ibd_m2hvl_gsego <- m2hvl_deg |>
  filter(cluster == 'IBD', p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', eps = 0)

ibd_m2hvl_gsego <- ibd_m2hvl_gsego |>
  simplify()

ibd_m2hvl_gsego@result |>
  filter(NES > 0, ONTOLOGY == 'BP') |>
  plot_enrichment(metric = NES) +
  labs(title = 'GO BP pathway GSEA in IBD colon:\nTRPM2-hi MF vs TRPM2-lo MF')

ibd_m2hvl_gsego@result |>
  write_source_csv('ibd.colon.mf.m2hvl.go.gsea')
